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ABSTRACT 

We present fully three-dimensional local simulations of compressible MRI turbulence 
with the object of studying and elucidating the excitation of the non-axisymmetric 
spiral density waves that are observed to always be present in such simulations. They 
are potentially important for affecting protoplanetary migration through the action 
of associated stochastic gravitational forces and producing residual transport in MHD 
inactive regions through which they may propagate. The simulations we perform are 
with zero net flux and produce mean activity levels corresponding to the Shakura 
& Sunyaev a ~ 5 x 10^'^, being at the lower end of the range usually considered in 
accretion disc modelling. We reveal the nature of the mechanism responsible for the 
excitation of these waves by determining the time dependent evolution of the Fourier 
transforms of the participating state variables. The dominant waves are found to have 
no vertical structure and to be excited during periodically repeating swings in which 
they change from leading to trailing. The initial phase of the evolution of such a swing 
is found to be in excellent agreement with that expected from the WKBJ theory 
developed in a preceding paper by Heinemann & Papaloizou. However, shortly after 
the attainment of the expected maximum wave amplitude, the waves begin to be 
damped on account of the formation of weak shocks. As expected from the theory the 
waves are seen to shorten in radial wavelength as they propagate. This feature enables 
nonlinear dissipation to continue in spite of amplitude decrease. As a consequence the 
waves are almost always seen to be in the non linear regime. 

We demonstrate that the important source terms causing excitation of the waves 
are related to a quantity that reduces to the potential vorticity for small perturba- 
tions from the background state with no vertical dependence. We find that the root 
mean square density fluctuations associated with the waves are positively correlated 
with both this quantity and the general level of hydromagnetic turbulence. The mean 
angular momentum transport associated with spiral density waves generated in our 
simulations is estimated to be a signiflcant fraction of that associated with the turbu- 
lent Reynolds stress. 
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1 INTRODUCTION 

This is the second of two papers which study the excita- 
tion of spiral density (SD) waves in accretion discs occur- 
ring through the action of stochastic forcing due to tur- 
bulence. In practice this has been taken to arise from the 
magneto-rotational instability (MRI) but excitation result- 
ing from stochastic forcing produced by other mechanisms 
is expected to lead to similar results. In the preceding pa- 
per ( [Heinemann fc Pap aloizou 2009j paper I) we developed 
a WKBJ theory of the excitation process. In this paper we 



study the wave excitation process directly as it is manifested 
in fully nonlinear three dimensional numerical simulations, 
with the object of elucidating it further and comparing the 
results with the predictions of the WKBJ theory developed 
in paper I. 

The outline of this paper is as follows. In Section [2] 
we give a brief review of the shearing box model, the ba- 
sic equations solved in the simulation, as well as the equa- 
tions governing the dynamics of SD waves. We also describe 
the procedure for performing the Fourier transforms of the 
state variables participating in the excited waves in shearing 
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coordinates. This allows us to identify and then follow the 
evolution of wave amplitudes in order to make a detailed 
comparison with the WKBJ theory presented in paper I. 

In Section [3j after briefly describing the simulation set 
up and the numerical method used, we present a fully three- 
dimensional simulation of compressible MRI turbulence in 
which we directly observe the excitation of large scale, non- 
axisymmetric SD waves. We discuss their main character- 
istics and the correlation of the root mean square density 
fluctuations associated with them with the level of turbu- 
lent activity. We go on to study the source terms driving 
the wave excitation and motivate the simplifying assump- 
tion that those proportional to a quantity that is related to 
the potential vorticity dominate. This assumption enabled 
us to derive the linear theory of wave excitation presented in 
paper I. A detailed comparison between this linear theory 
and simulation outcomes is presented in Section |4] There 
we follow the time dependent evolution of the appropriate 
Fourier amplitudes through a swing cycle, during which the 
wave excitation occurs, in detail. We obtain very good agree- 
ment with predictions from the theory developed in paper I. 
Finally, we discuss and summarize our results in Section [5j 
indicating the magnitude and the scaling of the wave angu- 
lar momentum flux with the Reynolds stress and also the 
wavelength in the azimuthal direction for which the wave 
excitation phenomenon is optimal. 



2 PRELIMINARIES 

2.1 The shearing box model 

As in paper I we consider an isothermal conducting gas in 
the shearing box approximation in the absence of vertical 
stratification. The Cartesian coordinate system (a:, y, z) is 
such that X corresponds to radius, y to azimuth, and z de- 
notes the vertical direction. The equations of motions consist 
of the continuity equation 

Vp + V-p^O, (1) 

the momentum equation 

Vp = -cVp -2flxp + qnp^By + V ■ T + V ■ {2puS) (2) 

and the induction equation 

VB = V X (u X B - 77V X B) - qflB^By (3) 

where the differential operator 

V = dt^qnxdy, (4) 

accounts for advection by the linear shear —q^lxBy, fl = fle^ 
is the angular velocity, p is the density, the linear momentum 
density p = pu defined in terms of the velocity deviations 
u = [ux, Uy, Uz) from the mean shear, c = HVl is the isother- 
mal sound speed, H is the nominal density scale height, the 
magnetic field B = {Bx, By,Bz) is normalised by the square 
root of the vacuum permeability pa, in which case the non- 
linear stress tensor T has components 



Tij = BiBj — SijB /2 — puiUj 



(5) 



S is the traceless rate-of-strain tensor whose components are 
given by 



5. . - 1 ( _|_ 

2 \dxj dxi 



which are seen to include terms arising due to the linear 
shear. The kinematic viscosity is v and the resistivity is rj. 

We consider the governing equations ([l]|3| to be subject 
to 'shearing periodic' periodic boundary conditions, i.e. 

f{x + Lx,y~ qntLx, z, t) = f{x, y, z, t), (6a) 
f{x,y + Ly,z,t) = f(x,y,z,t) (6b) 



and 



f(x,y,z + L^,t) = f{x,y,z,t). 



(6c) 



From Q we see that the domain becomes fully periodic once 
every shearing time 



ST, = 



Ly I Lx 

qO. 



(7) 



2.2 Wave equations 

We developed equations governing the excitation of SD 
waves in the inviscid limit in paper I. These follow directly 
from equations ([T| and ([2|. As we will demonstrate, SD 
waves in the unstratified shearing box exhibit little depen- 
dence on the vertical coordinate z so that we may consider 
vertically averaged quantities (indicated by being contained 
within angle brackets {■)z) for which the wave equations read 



(D^ _ + t^) (p)^ + 2qndy (px), = 

-2n{c)^^dx{mx),-dy{my)^, (Sa) 

p2 _ c2v" + k') {px), + 2qQc^dy (p)^ = 

c^dy{0^ + 2n{my)^+v{mx)^, (sb) 

and 

- c'dx (C), + {q- 2)n {mx)^ + V {^y)^ . (Sc) 

Here, = 2(2 - g)^^ is the square of the epicyclic fre- 
quency and we have introduced the divergence of the non- 
linear stress tensor = V ■ T as a short-hand to denote the 
nonlinear source terms. 

In paper I we postulated that as far as wave excitation 
is concerned the dominant source term in each of the three 
wave equations is the linear one appearing first on the right 
hand sides of ([sl and involving the quantity 



(C). = dx (Py), - dy (Px)^ + {q~ 2)n (p)^ 



(9) 



which we refer to as pseudo potential vorticity (or PPV for 
short). An important property of PPV is that apart from ad- 
vection by the linear shear, it varies in time due to nonlinear 
stresses only. 



v{o, = dx{^y),-dy{mx)^, 



(10) 



and is thus conserved for a linear SD wave. Note also that 
in this limit the change in PPV is equal to the change of 
potential vorticity, the difference being at most second order 
in the wave amplitude. 
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2.3 Decomposition into shearing waves 

In the periodic shearing sheet a suitable Fourier basis 
is given by plane waves with, in the coordinate system 
adopted, time dependent radial wave numbers. The expan- 
sion of any (vertically averaged) fluid variable / in this basis 
reads 

= f(t) exp [ik^{t)x + ikyy], (If) 



where / is a complex valued wave amplitude and where the 
wave numbers 



and k 

L/y 



Each of these so-called shearing waves ( |Thomson|p87| IS 
uniquely labelled by the pair of integral numbers n^jUy £ Z. 
The radial wave numbers depend linearly on time such that 



lit 



= qQky 



which is a consequence of advection by the linear background 
shear. 

Because kx{t)/ky increases monotonically if qi} > 0, 
every shearing wave evolves from being leading, i.e. 
kx{t) /ky < 0, to being trailing, i.e. kx{t)/ky > 0, as time pro- 
gresses from t = — oo to t = oo. The change from leading to 
trailing is referred to as 'swing' and occurs at 



ts = 



(12) 



From linear theory we expect wave excitation to occur at 



the time of the swing. From ( 12 I we see that different shear 



ing waves swing from leading to trailing at different times 
so that the excitation process consists of a series of swings 
that are separated, for a given ky, by the fixed time interval 

Sts = STs/Uy. 

Expanded in the Fourier basis ( |11[ ), the SD wave equa- 
tions ([8| for a single shearing wave read 



d^ 
dt2 



-I- [fe {t)c + k'^] p + 2qQikyPx 



2Q,C-ikxit)mx-ikymy, (13a) 



dt2 

and 
dt2 



-I- [k^{t)c^ + K^] px + 2qQc^ikyp = 



c"ifc„C + 2J79tj; + (13b) 



<mx 

dt 



+ [k'^(t)C^ + K^] Py = 



dTly 



cikxC + {q-2)nmx + -^. (13c) 



3 SIMULATIONS 

The steady state solution for the non-stratified shearing box, 
for which the density is uniform and the velocity takes the 
from of a linear shear, is subject to the magneto-rotational 
instability when a relatively weak magnetic field with zero 



net flux is introduced. This ultimately results in the pro- 
duction of turbulence that is powered by the linear shear 



and attains a statistically steady state (Hawley et al.|1995 



[Brandenburg et al.|1995| [Stone et al.|1996| . 

We have carried out local, three-dimensional simula- 
tions of accretion disk turbulence induced by the MRI. The 
purpose of these simulations is threefold. Firstly to study 
and quantify the mode of excitation of SD waves in a turbu- 
lent shearing box, secondly to justify some of the simplifying 
assumptions that were made in order to proceed with the 
theoretical analysis of the excitation mechanism given in 
paper I, and thirdly to provide results which can be tested 
against predictions made from this theory (see Section [4| so 
that the nature of the phenomenon and the implied depen- 
dence on physical parameters can be conflrmed. 



3.1 Setup 



We use the Pencil Code (see e.g. [Brandenburg fc Dobler[ 



2002 1 to solve the MHD equations in the shearing box. The 



numerical method relies on using sixth-order central finite 
differences for the evaluation of spatial derivatives and a 
third-order Runge-Kutta scheme for time integration. The 
method is stabilized by imposing explicit diffusion coeffi- 
cients in all evolution equations. The shearing box boundary 
conditions are implemented via sixth-order polynomial in- 
terpolation. The code solves the isothermal MIfD equations 
in the shearing box approximation ([l]-[3| in terms of mass 
density p, fluid velocity u, and, to ensure that V ■ B = 0, 
the induction equation is rewritten in terms of the magnetic 
vector potential A with B — V x A. Details are given in 
Appendix [A] 

The setup is similar to the one used for a code 
comparison in [Fromang et al.[ ( |2007| |. We choose a box 
size of {Lx,Ly,Lz) = {H,AH,H) and a resolution of 
{Nx,Ny,N^) = (128,512,128). [Fromang ( |2007 l 
showed that in the zero-net-flux case, the level of MRI tur- 
bulence critically depends on the magnetic Prandtl number 
Pm = v/r], i.e. the ratio of kinematic viscosity to magnetic 
resistivity. Here we set v = Z2-1Q~^ H^Q. and 77 = 8-10"^-ff^n 
such that Pm = 4. In this case we observe MRI turbulence 
that is, as in Fromang et al. (20071, sustained for several 



hundred orbits and expected to be numerically resolved for 
the values of the diffusion coefficients that we impose. 



3.2 Volume averages 

In Fig. [1] we plot the root mean squared velocity and density 
fluctuations, calculated as volume averages over the box, 
and volume averages of the Reynolds and Maxwell stress as 
functions of time. The initial peak that as seen in the velocity 
fluctuations and in the Reynolds and Maxwell stresses is due 
to exponential growth of axisymmetric MRI modes. After 
roughly 4 Orbits, this so-called channel solution is torn apart 



by secondary, non-axisymmetric instabilities (see e.g. Latter 
[et al. 2009, and references therein) and the system enters a 
quasi steady state of sustained MRI turbulence. 

The velocity fluctuations in this state are an order of 
magnitude less in amplitude than the speed of sound. In 
spite of the strongly subsonic nature of the turbulence we 
observe density fluctuations with a temporal mean of about 
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Figure 1. Temporal evolution of root mean squared velocity and 
density fluctuations (upper two panels) and of the volume aver- 
aged Reynolds and Maxwell stress (lower two panels). The red 
dashed line corresponds to t = 82 5Tb at which point the largest 
density fluctuation occurs. 



6 per cent and peak values of more than 15 per cent of the 
background value. We will see that the peak density fluctua- 
tions are due to high-frequency, non-axisymmetric SD waves 
that are excited via the mechanism discussed in paper I. Be- 
cause these waves are trailing, they lead to outward angu- 
lar momentum transport. This manifests itself in a strong 
correlation between the density fluctuations and the mean 
Reynolds stress. At peak values, the Reynolds stress may 
even become comparable to the Maxwell stress although we 
note that the former fluctuates rapidly in time so that the 
net angular momentum transport due to SD waves is less sig- 
nificant than it may appear at first glance. We will attempt 
to quantify the amount of angular momentum transport due 
to SD waves later in this paper. 

We note that the Maxwell stress fiuctuates in time much 
less rapidly than, and is only very weakly correlated with, 
the other (purely hydrodynamic) quantities. This is because 
SD waves are non-magnetic and the Maxwell stress varies 
on a time scale comparable to the turn over time scale of 
the turbulence which in this case is much longer than the 
inverse SD wave frequency. 
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Figure 2. Pseudo-colour image of mass density on the faces of 
the computational domain. This snapshot is tal^en at t = 82 STs 
where we observe the largest peak in the volume averaged density 
fluctuations shown in Fig. [l] 



3.3 Flow structure in real space 

The typical large scale three-dimensional structure of the 
turbulent fiow is depicted in Fig. [2] where we show the mass 
density field projected onto the faces of the computational 
domaiij^ This snapshot is taken at f = 82 STs where we ob- 
serve the largest peak in the volume averaged density fiuc- 
tuations shown in Fig. [l] 

The wave crests are seen to propagate at an inclined 
angle within the horizontal xy-plane and display little de- 
pendence on z. In order to ease the subsequent analysis, as 
in paper I, we may thus integrate over z and will look at 
vertically averaged quantities only from now on. 

In Fig. [3] we plot false colour images of such vertical av- 
erages at the same instance in time as shown in Fig. [2] Large 
scale wave crests are clearly visible in (p)^, {Px}^ and, to a 
lesser extent, in (Py)^- The waves are trailing and always 
occur in pairs where one wave propagates in a direction op- 
posite to the other wave. That this has to be the case is clear 
from symmetry considerations of the shearing sheet and we 
recall that our linear theory predicts pairwise excitation of 
waves with equal amplitude but oppositely directed wave 
vectors. 

In contrast to (p)^, {Px)^, and (py)^, the pseudo poten- 
tial vorticity 

(C), = dx (Py)^ - dy (px), + (g - 2)n {p}^ 

does not display a corresponding large scale wave structure 
but instead consists of smaller scale turbulent structures 
that are, as is characteristic for accretion disc turbulence, 
elongated in the shearwise direction. The absence of large 
scale wave structure in {(^) ^ is reminiscent of the fact that 
{C)^ is conserved for a linear wave, see (10 1, and the justifi- 
cation for why we are allowed to treat {q7 as a source term 
for wave excitation, as was done in paper I, in the first place. 



^ An animated version of Fig. [2] and Fig. [3] below is provided as 
part of the on-line supplements to this article 
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Figure 3. Pseudo colour images of z-averaged quantities at the same instance in time as in Fig. |2] Non-axisymmetric waves are clearly 
visible in {p) ^ and {px) and to a lesser extent in {py) ^- There will be a system of regularly repeating wave crests when boxes are stacked 
together to fill space by application of the periodicity conditions. In contrast to this, the waves are not visible in {&Cl where SC, = C, — C,q 
is the PPV deviation from the constant background value Co = (</ ~ 2)f!po. 



3.4 Evolution of single shearing wave 

A more revealing picture of the evolution of SD waves as 
observed during the simulation may be obtained in Fourier 
space. For this purpose we will look at the discrete Fourier 
transform (DFT) of vertically averaged quantities as shown 
in Fig. [3] Taking a two-dimensional spatial DFT in the 
shearing sheet is non-trivial due to the presence of the linear 
background shear and requires a coordinate transformation 
to shearing coordinates, see Appendix |B] This coordinate 
transformation can be done very conveniently by first tak- 
ing the DFT along the periodic j/-direction, then applying 
an x-dependent phase shift to achieve full periodicity in x, 
and finally taking the DFT along the a;-direction. Further 
details are given in the Appendix |B] This procedure allows 
us to uniquely assign a time-dependent radial wave number 
kx(t) as well as a (constant) azimuthal wave number ky to a 
given Fourier amplitude obtained from the DFT, enabling us 
to follow individual plane waves as they swing from leading 
to trailing during the course of the simulation. 

In Fig.|4]we plot the evolution of the wave amplitudes of 
the p, Px, and py for a single shearing wave with azimuthal 
wave number ky = 2ti/ Ly. Excitation occurs when the wave 
swings from leading to trailing, i.e. when kx{t) « 0, which for 
this particular wave happens aX t = S15Tb. The maximum 
density amplitude of Re p ~ 0.08 po is obtained when kx (t) ~ 
2n/H or, equivalently, around t = 82STs, i.e. roughly one 
shearing time later. In Fig. [4] we are thus looking at the 
Fourier representation of the wave seen in Fig. [2] and Fig. [3] 
and thus at that particular wave responsible for the large 
peak in the volume averaged density fluctuations indicated 
by the dashed red line in Fig. [T] 



We note that the wave undergoes significant damping 
for kx{t)H > 2n. This decay is in stark contrast to linear 
theory which either predicts an algebraic decrease (for py) 
or even increase (for /} and px) of the wave amplitude with 
time (see paper I) and therefore needs to be attributed to 
nonlinear effects. We will comment on this further in Sec- 
tion H 

The simulation allows us not only to monitor the wave 
amplitudes p, px, and Py for a single shearing wave as it 
evolves in time, but also the various source terms occurring 



on the right sides of the SD wave equations (131, and to 



assess the relative importance of these terms for wave exci- 
tation. 

There are both linear source terms involving the pseudo 
potential vorticity as well as nonlinear terms involving the 
vector Vt = ik T, where Tij is the Fourier transform of the 
nonlinear stress tensor (|5|. 

In Fig. [5] we plot the Fourier amplitudes of ^ and the Tij 
for the wave shown in Fig. [4] We see that C is - compared to 
the p, Px, Py, and the nonlinear stresses Tij - a slowly vary- 
ing function of time. When wave excitation occurs around 
kx{t) = there is no sign of oscillating behaviour in This 
is a manifestation of PPV conservation to linear order, see 

ITol. 



Furthermore, at the time of wave excitation, the ampli- 
tude of is more than an order of magnitude larger than the 
amplitude of any of the components of the nonlinear stress 
tensor Tj suggesting, as these quantities measure their rel- 



ative contributions to the source terms in equations (131, 



that it is indeed the pseudo potential vorticity, and thus the 
linear source term, that is primarily responsible for wave 
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Figure 4. Fourier amplitudes of density and linear momenta for a 
single shearing wave with azimuthal wave number ky = 2-k/ Ly as 
it swings from leading to trailing. Excitation oecurs when kx (t) = 
or, equivalently, at t = 81 ST^. This is thus the same wave that 
is shown in Fig. [2] and Fig.[3] 



Figure 5. Fourier amplitudes of the pseudo potential vorticity 
and the components of the nonlinear stress tensor. All of these 
appear as sources on the right hand side of the SD wave-equations 



excitation. We will provide more evidence for this claim in 
Section ID 



4 COMPARISONS WITH ASYMPTOTIC 
THEORY 

The evolution of SD waves in the shearing sheet is de- 
scribed mathematically by the time-dependent forced oscil- 



lator equations (13 1. Wave excitation is due to the presence 



of inhomogeneous source terms in these equations and is 
found - as demonstrated in the preceding section - to occur 
at the swing stage, i.e. when the time dependent radial wave 
number kx{t) = 0. 

In paper I we derived analytic solutions to the governing 



equations ( 13 1 that are able to capture this swing excitation. 



These solutions are asymptotic in the parameter 

qQkyC 



which is small in both the low and the high azimuthal wave 
number limit, given by fcyC ^ k and kyC^ k, respectively. 
The asymptotic solutions can be written as the sum of a non- 
oscillatory, vortical solution - called the balanced solution - 
to the inhomogeneous problem, and standard WKBJ solu- 
tions to the free wave equations that contain the oscillatory 
part of the solutions and describe the excited SD waves. 

A key hypothesis in our analysis in paper I was that the 
dominant source term for wave excitation is the linear term 
involving the pseudo potential vorticity C, and that nonlin- 
ear stresses contribute to the excitation only indirectly by 
generating C, prior to the swing via (10 1. This assumption 



kyC^ + K? ' 



renders the analysis as essentially linear. We have shown in 
Section |3.4| that the pseudo potential vorticity for a large 
amplitude SD wave observed in the simulation is more than 
an order of magnitude greater in size than the nonlinear 
stresses when the excitation occurs. This suggests that our 
quasi-linear description of the excitation process is indeed 
valid and we are now in a position to verify this assertion 
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by comparing the observed wave amplitudes with our pre- 
dictions from hnear theory. 

Expressed in dimensional form, the leading order bal- 
anced solutions for p, px, and py can be obtained from the 
relations given in paper I as 



Px = iCcIm 
Py — i^c Re 



ikyC — 



2qQ,\kyC 
- 2n 

k^c^ + + 2qD,ikyC 
ikxc 



\kyC - 



where the SD wave frequencies 



(14a) 
(14b) 
(14c) 



and 



Lip- ^ 2qSX\kyC. 



Note that these solutions include the full time dependence 
of the source. 

We are, however, mainly interested in the oscillatory 
wave part of the solutions, i.e. the WKBJ solutions to the 
homogeneous problem given in paper I, which may be writ- 
ten as 



-C,Sl ^Im 



LU^ dt' 



and 



Px = iCs-ffRe 



oj-i- sm 



/ loj^ dt' 
Jo 



(15a) 



(15b) 



where 



A, 



(A:2c2 + «:2 + 2qnikycy/* 



2n 



ikyC 



^■/4e 



(16) 



The free WKBJ solution for py follows from PPV conserva- 
tion, 



Py = 



ikyPx + 2(2 — q)Q,p 
\kx 



In the above expressions, only the value of the pseudo po- 
tential vorticity at the time of the swing (denoted by C,b) 
enters, to which the amplitude of the excited wave is di- 
rectly proportional. This amplitude is exponentially small 
in the asymptotic limit e ^ 1. However, we demonstrated 
in paper I that our leading order asymptotic solutions are 
remarkably accurate even for e < 1, and in the case of Kep- 
lerian shear effectively hold for the entire range of azimuthal 
wave lengths because e ^ 3/4 for ^ fc^i? < oo. 

4.1 Comparison of wave amplitudes 

We are now in a position to test the linear theory of the 
excitation process by comparing the asymptotic solutions 
\\A\ and ( |15[ ) with the evolution of an individual shear- 
ing wave observed in the simulation as shown in Fig. |4] and 
Fig. [5] For this particular wave, the azimuthal wave number 
ky — 2tx /Ly and thus e « 0.68. 

The result is shown in Fig. [6] As in the linear problem 
discussed in paper I, the numerical solution in the leading 
phase follows the balanced solution very closely. In the trail- 
ing phase the WKBJ solution provides an excellent match 




Figure 6. Comparison between simulation and linear theory for 
tlie same shearing wave as in Fig. |4] The curves shown in each 
panel are the wave amplitude as recorded during the simulation 
(black), the balanced solution (blue), and the full WKBJ solu- 
tion (red). For each variable, the WKBJ amplitude in the upper 
panel is taken as is from linear theory and therefore grows/decays 
algebraically with time. 



to the numerical solution for small radial wave numbers 
kx < 2-k/H. Roughly one shearing time after excitation, 
when the radial wave length is equal to the scale height H , 
the numerical solution starts to undergo significant damp- 
ing whereas the amplitude of the WKBJ solution continues 
to increase algebraically. Empirically we find that the de- 
cay of the numerical solution cannot simply be attributed 
to viscous effects because the diffusion coefficients imposed 
for this simulation are at least an order of magnitude too 
small. 

We are thus lead to believe that SD waves are damped 
due to nonlinear effects. It is clear that a linear shearing 
wave whose radial wave length decreases as its amplitude 
increases eventually has to break and form a shock. This 
phenomena has been already discussed in the context of disk 
physics applied to e.g. Saturn's rings ('Goldreich & Tremaine 
|T978| ) and to spiral density wakes of protoplanets ( Goodman 
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Figure 7. Upper panel: Volume averaged density fluctuations as 
a function of time. Here we have used high order binomial smooth- 
ing to get rid of the fast oscillations seen in Fig. [l] Lower panel: 
Modulus of the potential vorticity amplitude for the two smallest 
non-zero azimuthal wave number in the system {ky = 2-n jLy and 
ky = -k/ Ly). Each red dot corresponds to an instant in time when 
there is a shearing wave with k^it) = 0. This happens once (twice) 
per shearing time for waves with ky = 1-k jLy [ky = ir/Ly). 



\&i Rafikov||200l] ). If damping is mostly due to this efTect, 
we conclude from Fig. [6] that shock formation must happen 
already when the radial wave length is roughly equal to H. 
This fits with the observation that SD waves usually appear 
as weak shocks in the simulation, see Fig. [2] and Fig. |3] as 
well as e.g. Gardiner & Stone (20051. 



4.2 The importance of PPV 

The results from the previous subsection indicate that lin- 
ear theory is able to account for the initial evolution of an 
excited SD wave observed in the simulation. From ( |15[ ) we 
see that the only dynamic quantity that enters the linear 
SD wave amplitude is the pseudo potential vorticity at the 
swing state, Ca- K the particular wave discussed in the previ- 
ous subsection was not just a lucky coincidence, and linear 
theory holds in general as far as the initial evolution is con- 
cerned, we thus expect a correlation between (^s and, for 
instance, the root mean squared density fluctuations. 



shown in Fie. [T] 

In Fig.Mwe plot ((Sp)rms as well as j(Js| for the two short- 
est, non-zero azimuthal wave numbers available in our com- 
putational domain, denoted here by ky and fc", as a function 



of time. The temporal correlation between (5p)ims and |(^s| 
for ky — fcy is remarkably strong indeed. In particular the 
peak density fluctuations are seen to be delayed by roughly 
one shearing time with respect to the corresponding peak 
in |(^s|. This is because it takes roughly one shearing time 
for an SD wave to amplify to full strength after excitation, 
see Fig. [6] This hardly obvious correlation is perhaps the 
strongest evidence for the validity of our linear theory and 
the neglect of nonlinear source terms in the calculation of 
the wave amplitude. 

In contrast to the above, the temporal correlation be- 
tween (5p)rms and |<Js| for ky = k^ is much weaker even 
though the average magnitude of is still roughly 70 per 
cent of the average magnitude for ky ~ k^. On one hand, this 
is likely to be the case because the correlation for ky — k^ is 
obscured and therefore harder to see, but more importantly 
because swing excitation - as we know from linear theory - 
is simply a weaker effect for ky = k]^ than for ky = ky due to 
the exponential dependence of the WKBJ amplitudes on ky, 
see (16 1. We will comment on this further in the discussion 



section. 



4.3 Angular momentum flux 



In paper I we derived an expression for the radial angular 
momentum flux associated with a single pair of linear SD 
waves propagating in opposite directions. 



-2kyC Im(^*p 



(17) 



see Goodman & Ryu. The angular brackets denote an av- 
erage over azimuth, and as in paper I we take ky 0. The 
radial Lagrangian displacement ^a, is defined via 



dt 



Po' 



Note that (171 only includes the oscillatory wave part of the 



Fourier amplitudes, i.e. we do not account for the supposedly 
small contributions to the total angular momentum trans- 
port resulting from the slowly varying balanced solutions 
fT4|. 

Since from a numerical point of view it is difficult to 
unambiguously separate the oscillatory wave part of ^x, it 
is - for diagnostic purposes - more convenient to use the 
modified angular flux 



y- 2(2 - q)Qpo 



Im(p*p) 



(18) 



which does not involve the Lagrangian displacement. The 
two alternative expressions for the angular momentum flux. 



( 17 1 and ( 18 1, are found to agree when averaged in time over 
one oscillation period. 

It is now of interest to compare the angular momentum 



flux ( 18 1 associated with a single pair of SD waves observed 



in the simulation with linear theory by substituting the free 



wave solutions (151 into (181. The wave part of the Fourier 



amplitudes obtained from the simulations is calculated by 
subtracting the balanced solutions ( |14[ ). 

The result is shown in in Fig. |8| where we plot the mod- 
ifled angular momentum flux associated with the SD wave 
displayed in Fig. [6] As expected there is essentially no an- 
gular momentum transport prior to excitation, i.e. as long 
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Figure 8. Angular momentum flux associated with the SD wave 
displayed in Fig. |6] 



the waves are leading. In the trailing phase, after swing ex- 
citation has occurred, the angular momentum flux oscillates 
rapidly due to interference between the forward and the 
backward travelling wave. However, it does so around a pos- 
itive mean, which means that there is net outward angular 
momentum transport. 

In contrast to linear theory, the angular momentum flux 
observed in the simulation is strongly damped for moder- 
ately large wave numbers due to nonlinear effects. Interest- 
ingly enough, however, before the onset of damping the flux 
for a single pair of shearing waves observed in the simulation 
reaches a maximum value that is very close to what is pre- 
dicted by linear theory for the far trailing regime. In spite 
of nonlinear damping, there will be a steady flux of angular 
momentum due to successive swings of shearing waves and 
we can thus use linear theory to make an order of magni- 
tude estimate of the angular momentum transport due to 
SD waves. 

In linear theory we calculated the net angular momen- 
tum flux in the far trailing regime explicitly in paper I. Aver- 
aged over one oscillation period, this flux attains a constant 
non-zero value given by 



[(fc2c2 + k2)2 + {2qVlkyC) 



i3/4 ' 



(19) 



as kxit) oo. In the above equation. 



tan"^(2e) 



2e 



In order to estimate the amount of angular momentum 
transport due to SD waves in the simulation, we thus re- 
quire the swing PPV amplitude ^g. For shearing waves with 
the longest azimuthal wave length available in our computa- 
tional domain, this information can be extracted from Fig. [7] 
where we plot the modulus of the swing PPV amplitude for 
all shearing waves that we observe to swing from leading 
to trailing. Making an ensemble average over all swings dis- 
played in Fig. [7|we obtain 



lo^^n^Po- 



Substituting this into (|19| yields 



lim {Fx) 

kx/ky — >-oo '^^ 



2.5 ■ 10~ poc 



From this rough estimate we conclude that the mean an- 
gular momentum transport due to SD waves in the present 
simulation is possibly less but still significant compared to 
the turbulent Reynolds stress, see Fig. [l] 



5 SUMMARY AND DISCUSSION 

In this paper we have presented fully three-dimensional sim- 
ulations of compressible MRI turbulence in which we ob- 
served the excitation of large scale, non-axisymmetric SD 
waves. The root mean square density fluctuations associ- 
ated with these waves were shown to be correlated with 
the degree of turbulent activity measured by the [Shakura] 

6 Syunyaev ( 1973 1 a parameter. 



(Bx By 



pUxUy)^ 



PQC^ 

We performed a detailed analysis of the excitation of 
these waves by following the time dependent evolution of 
the Fourier transforms of the participating state variables. 
Their evolution during the excitation phase of the associ- 
ated SD wave was found to be in excellent agreement with 
the WKBJ theory developed in paper I. However, shortly 
after the attainment of the maximum wave amplitude, the 
waves are rapidly damped. This damping is consistent with 
the sharp density profiles observed in the simulations which 
are indicative of the existence of shocks. Shocks are expected 
even under conditions of weak excitation because as shearing 
waves propagate their radial wavelength shortens, resulting 



in the onset of nonlinearity and shock formation (e.g. Good- 
man & Rafikov 2001 1. Thus the waves are always likely to be 



observed in the nonlinear regime. Our results indicate that 
for activity levels corresponding to q ~ 5 x 10~^, being at 
the lower range of magnetic activities normally considered, 
shock formation occurs when their radial wavelength is com- 
parable to the scale height. 

The excited waves are trailing and associated with 
an outward angular momentum flux. The calculations per- 
formed here are consistent with the WKBJ theory presented 
in paper I in indicating that this outward flux scales in the 
same way as the Reynolds stress. This is supported by the 
fact that the density fluctuations associated with the waves 
are correlated with it. The level of angular momentum trans- 
port associated with SD waves is estimated to be a signif- 
icant fraction of that due to the turbulent Reynolds stress 

{~pUxUy)y. 

By consideration of the source terms for the wave exci- 
tation, we identified those associated with the pseudo poten- 
tial vorticity (PPV) as the dominant ones as was assumed 
in paper I. Wave like structure was not seen in the PPV 
field, see Fig. [Sj indicating that this could act consistently 
as an independent source generated by background turbu- 
lent fluctuations. 

For the simulation discussed in detail, SD wave excita- 
tion is found to be most effective for waves with the smallest 
azimuthal wave number available in the computational do- 
main. This can be understood from the theory presented 
in paper I. There it was found that the amplitude of an 
excited SD wave is directly proportional to the Fourier am- 
plitude of the pseudo potential vorticity at the swing stage, 
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vertical dependence. Although we considered an unstratified 
disc, SD waves with weak vertical dependence also exist in 



1/4 1/2 



Figure 9. Tiie spectrum of swing pseudo potential vorticity av- 
eraged over all swings as a function of azimuthal wave number ky . 
The spectrum peaks at the smallest ky, remains rather shallow 
at moderately small ky, but falls off for high ky. 



("s, while also having an exponential dependence on the az- 
imuthal wave number ky through the parameter 

qQkyC 



If the latter dependence alone is important it implies that, 
for a Keplerian rotation law, SD wave production will be 
most effective for ky ~ k"^^ — fl/c. This is longer than the 
longest possible azimuthal wavelength in the shearing box 
considered here with Ly — AH, so that SD wave production 
should most effective at the longest azimuthal wavelength. 

The other important dependence is that of the PPV 
spectrum on ky at the time of swing. This spectrum is diffi- 
cult to predict theoretically, but easily determined from the 
simulation data. We do so by taking for each ky an average 
over the ensemble of shearing waves that we find to swing 
from leading to trailing during the course of the simulation. 
The result is shown in Fig. |9] We see that the spectrum is 
rather shallow on large scales and falls off at small scales. A 
nearly flat spectrum at small ky means that the dependence 
of the wave amplitude on ky, at the most effective values, 
is determined by the WKBJ exponential factor appearing 
in e.g. ( |19[ l, and should be maximal at either ky — 1/H, or 
the smallest possible ky in the box should that be larger. 
Accordingly we expect the smallest ky to be dominant for 
Ly = AH but this dominance should begin to be lost once 
Ly starts to exceed 2-nH. We checked that when Ly was 
extended to 8//, the contribution from the smallest ky no 
longer exceeded that from the next smallest confirming the 
discussion given above. 

Differentially rotating discs undergoing MRI driven tur- 
bulence are ubiquitous in astrophysics and the results pre- 
sented here and in paper I indicate that the production of 
density waves is generic. The density perturbations associate 
with these may play an important role in driving protoplane- 
tary migration (e.g. [Nelson fc Papaloizou|2004[|Nelson|2005[ ) 
and may play a role in exciting quasi periodic oscillations in 
discs around close binary stars. In the case of protoplane- 
tary discs it has been suggested that only the high altitude 
regions may be turbulent while the mid-plane region lies in a 
'dead zone' (e.g. |Gammie|1996l ). We comment that we found 
that the SD waves excited by the turbulence have very weak 



the isothermal stratified case (e.g. Fromang & Papaloizou 
|2007[ ), so that although they may be excited in the upper 
layers propagation throughout the vertical structure is to be 
expected with the consequence that some small but non zero 
transport may result throughout the disc. In this context we 
note that, consistently with these ideas, recent simulations 



by Oishi et al. (20071 that incorporated dead zones were 
found to have wavelike motions excited in the mid-plane re- 
gions associated with a small residual angular momentum 
flux that was uniformly distributed in the dead zone with 
a typical a value between 10~^ and 10~^. We remark that 
a simple application of the WKBJ theory developed in pa- 
per I suggests, since the source of wave excitation is a ver- 
tical average of the pseudo potential vorticity, and the wave 
angular momentum flux is the square of this, that the value 
of ct in the dead zone should scale as ratio of the square of 
the surface density in the high altitude turbulent regions to 
the surface density of the whole disc. However, this may be 
a lower bound because the dependence on surface density 
may actually be weaker on account of weaker dissipation in 
the linear regime allowing the build up of a relatively larger 
amount of wave angular momentum flux. This aspect re- 
mains to be investigated and decided by future much higher 
resolution simulations that we are unable to perform at the 
present time. 
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APPENDIX A: THE PENCIL CODE 



The Pencil Code (see e.g. [Brandenburg fc Dobler|[2002l ) is 
a high-order finite-difference MHD code that is primarily 
designed to deal with weakly compressible turbulent flows. 
Solenoidality is ensured by solving the induction equation 
in terms of the magnetic vector potential. The equations of 
motion are evolved using an explicit Z^'^ order Runge-Kutta 
time stepping scheme. The scheme is stabilized by explicit 
diffusion coefficients. 

For the purposes of the present paper, the shearing box 
equations that the code solved numerically consist of the 
continuity equation 

T)p + u-Vp = — pV ■ u + pdiff, 

the momentum equation for an isothermal equation of state 

"Du + u ■ Vu — —2ft X u + qQuxBy 

(?Vp j X B 



UdiS, 



and the induction equation 

VA = qQAyEx + u X B + Adis, 
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where 



X. We can however transform to shearing coordinates, 



V — dt ~ qQxdy, 

p is the mass density, u is the fluid velocity, A is the mag- 
netic vector potential, B — V x A is the magnetic field, 
j = VV ■ A — \7^A is the current density, ft — Q,ez is the 
angular velocity, and c is the speed of sound. For Keplerian 
rotation q = 3/2. 

The numerical scheme is stabilised by imposing explicit 
diffusion terms. For the momentum equation and the induc- 
tion these take the same form as molecular diffusion, i.e. 



Wdiff = ( V w + - VV • u H 

3 P 



and 

4diff = V (V'A - VV • A) , 
where the traceless rate-of-strain tensor 



1 

dxj dxi 



-5jj V ■ u 



qQ. 



the values of kinematic viscosity v and magnetic resistivity 
are given in the text. 

In the case of the continuity equation, where there is no 
explicit physical diffusion coefficient, we maintain numerical 
stability by employing 6**^ order hyper diffusion, i.e. 



Pdiff 



where the mass diffusion coefficient S is chosen to scale with 
the grid spacing 5x, 

(Sxfc 



60 



so that significant mass diffusion only occurs near the grid 
scale. We note that the mass diffusion term ( [A| ) formally 
conserves mass and that we indeed find very little variation 
in the total mass during the course of the simulation. 

All derivative operators are given by 6*'' order accurate 
central finite-difference formulae, 



so that 



X — X 

y — y + qfltx 
dx = di + qQtdy 

dy = dy 

dt — df + qflxdy . 

In shearing coordinates ( |Bf [ | reads 

dif^c{ds + qntdy)f (B2) 

where the transformed field / is defined through 

/(i, y, i) = /(a:, y + q^tx, t) = f{x, y, t) (B3) 

and so is obtained by shifting y by qQ.tx. By transforming 
to shearing coordinates we have eliminated the explicit x- 
dependence at the expense of an explicit t-dependence on 



the right hand side of (B2|. We now impose fully periodic 



boundary conditions in shearing coordinates, i.e. 

/(£ + Lx,y..i) = f{x,y,i) 

f{x,y + Ly,i) = f{x,y,i). 
From ( |B3[ ) it then follows that 

f{x + Lx,y- qntLx , t) = f{x, y, t) 

and 

f[x,y + Ly,t) = f{x,y,t) 

which are just the shearing sheet boundary conditions that 
require periodicity in shearing coordinates 

We are now in a position to Fourier decompose any field 
f{x,y,t) in shearing coordinates. 



f{x,y,t)= ^ /n,,n„(t)exp 



2'Knx\ - . (l-KUy 
II ^Ix-I- 



where the Fourier coefficients can be computed from the 
inverse transform 



45(/.+l - - 9(/.+2 - .A-2) + (.A+3 - /»-3) 

60(52; fn^,ny{i) 

-490/, + 270(/,+i + - 27(/,+2 + ,^-2) + 2(/,+3 + .A-3) 



LxLy 



(6) _ 20/, - f5(/,+l + /,-l) + 6(/,+2 + f.-'z) - (/,+3 - f^-i) 



exp 



. ( 2-Knx 



dx dy. 



fr = 



{Sx) 



The shearing periodic boundary condition ( |6a[ ) is im- 
plemented using 6"^ order polynomial interpolation. 



APPENDIX B: DISCRETE FOURIER 
TRANSFORMS IN THE SHEARING SHEET 

Let us consider a prototypical evolution equation for a quan- 
tity f(x, y, t) in the shearing sheet of the form 



\ Lx 

Because dx dy = da; dy it follows from 

W ^ ^ //,/(.,., exp [- i.. M d, 

where the wave numbers in unsheared coordinates are 

' 2?™ 



, 2'n:nx \ , 2imy 

kx{t) = ( + i^^y-j;^ 



and ky = ( 

Ly 



More explicitly 



{dt - qflxdy) f ^ cdxf 



(Bl) 



We may Fourier-decompose this equation in y, but because 
of the explicit a:;-dependence we are not allowed to do so in 



exp 



f{x,y,t) exp 



—iqflt 
, / 2Tvny 



2Tmy 



dy } exp 



2TTnx 
Lx 



dx 



In order to compute the Fourier transform of a real 
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quantity from the simulation we therefore do an FFT in y 
first, muhiply the resuh with 



exp 



2-Kny 



(B4) 



and finally do an FFT in x. We then determine the time 
dependent radial wave number that each Fourier amplitude 
corresponds to from 



/ 2Tvnx 



+ qQ.t 



2Tmy 

Ly 



ttNx 



where Nx is the total number of grid points in the a;-direction 
and the integer m is chosen such that we always have 



T^Nx 

Lx 



<: kxit) <: 



ttNx 
Lx 



Note that the phase shift ( B4 1 corresponds to an x- 
dependent translation along y by an amount —qQtx which 
makes the field fully periodic in x. 
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